In [24]:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.''')
Out[24]:
The raw code for this IPython notebook is by default hidden for easier reading. To toggle on/off the raw code, click here.
In [ ]:
import rpy2
print(rpy2.__version__)
import rpy2.rinterface
from rpy2.robjects.packages import importr
import rpy2.robjects.packages as rpackages
#%load_ext rpy2.ipython

import rpy2
import rpy2.robjects.packages as rpackages
import pandas as pd

# import rpy2's package module
# import R's utility package
utils = rpackages.importr('utils',lib_loc="C:/Program Files/R/R-3.5.0/library")
#importr("ggplot2", lib_loc= Where_is_R)
#utils.chooseCRANmirror(ind=1) # select the first mirror in the list
#base=rpackages.importr('base')
# select a mirror for R packages

#import R's "base" package


!curl -O "https://raw.githubusercontent.com/vitorcurtis/RWinOut/master/RWinOut.py"
%load_ext RWinOut

###Python Multiplotter

def showImagesMatrix(list_of_files,hSize, wSize,col=10):
    fig = figure( figsize=(wSize, hSize))
    number_of_files = len(list_of_files)
    row = number_of_files/col
    if (number_of_files%col != 0):
        row += 1
    for i in range(number_of_files):
        a=fig.add_subplot(row,col,i+1)
        image = imread(mypath+'/'+list_of_files[i])
        imshow(image,cmap='Greys_r')
        axis('off')
        
from IPython.display import HTML

#HTML('''<script>
#code_show=true; 
#function code_toggle() {
# if (code_show){
 #$('div.input').hide();
 #} else {
 #$('div.input').show();
 #}
 #code_show = !code_show!curl -O "https://raw.githubusercontent.com/vitorcurtis/RWinOut/master/RWinOut.py"
%load_ext RWinOu
#} 
#$( document ).ready(code_toggle);
#</script>
#<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')

Southern Appalachian LANDIS-II parameter work

In [ ]:
!curl -O "https://raw.githubusercontent.com/vitorcurtis/RWinOut/master/RWinOut.py"
%load_ext RWinOut
In [17]:
def showImagesMatrix(list_of_files,hSize, wSize,col=10):
    fig = figure( figsize=(wSize, hSize))
    number_of_files = len(list_of_files)
    row = number_of_files/col
    if (number_of_files%col != 0):
        row += 1
    for i in range(number_of_files):
        a=fig.add_subplot(row,col,i+1)
        image = imread(mypath+'/'+list_of_files[i])
        imshow(image,cmap='Greys_r')
        axis('off')
In [16]:
from IPython.display import Image
PATH = "E:/FS_meeting/FiguresMisc/"
Image(filename = PATH + "Apps_Long.jpg", width=1000, height=1000)
Out[16]:
In [26]:
Where_is_R="C:/Users/zacha/Documents/R/win-library/3.5"
In [27]:
#import rpy2's package mod
importr("rgdal", lib_loc= Where_is_R)
importr("raster", lib_loc= Where_is_R)
importr("sf", lib_loc= Where_is_R)
importr("magrittr", lib_loc= Where_is_R)
Out[27]:
rpy2.robjects.packages.Package as a <module 'magrittr'>

C:N ratios and Lignin

Most of the nitrogen, carbon and lignin parameters were started from existing papers, and the TRY database. The file NECN_Folder contains the data taken from each paper and the database. Each Try is linked to an individual paper in the try files. Some species were adapted from existing papers where the formula was assumed to have a carbon ratio of 1/2. This formula is therefore C: N=.5 /(g nitrogen/g total). For species that were not included in either papers nor the plant database were researched individually, or given a value based on a genus/family/order similarity to a species in this or another Landis papers. I Started with values from previous LANDIS-II models and supplmented them with data from the TRY. Try data was queried by Species and by trait. In the even of multiple returns values were meaned.

TRY Data: Kattge, J., Bönisch, G., Günther, A., Wright, I., Zanne, A., Wirth, C., Reich, P.B. and the TRY Consortium (2012) TRY - Categorical Traits Dataset. Data from: TRY - a global database of plant traits. TRY File Archive https://www.try-db.org/TryWeb/Data.php#3

In [29]:
%%R
w_dir<-"E:/FS_meeting/NECN_Biochem/"
setwd(w_dir)
In [30]:
w_dir="E:/FS_meeting/NECN_Biochem/"
Try=pd.read_csv(w_dir+"Try_Proc.csv")
#Try.head(n=10)
print("Example Acer rubrum")
Try[Try.sp=="Acer rubrum"]
Example Acer rubrum
Out[30]:
Unnamed: 0 sp trait unit value
4 1 Acer rubrum Coarse woody debris (CWD) lignin content per C... NaN 77.750000
5 9 Acer rubrum Leaf nitrogen (N) content per leaf dry mass NaN NaN
6 16 Acer rubrum Branch coarse woody debris (CWD) nitrogen (N) ... % 0.264200
7 23 Acer rubrum Leaf nitrogen (N) content per leaf dry mass % 1.752372
8 59 Acer rubrum Stem (wood) carbon (C) content per stem dry mass % 48.640000
9 73 Acer rubrum Crown (canopy) width cm 2.812360
10 98 Acer rubrum Leaf carbon/nitrogen (C/N) ratio kg C per kg N 24.700000
11 111 Acer rubrum Leaf carbon (C) content per leaf area kg C per m2 39.900000
12 137 Acer rubrum Leaf nitrogen (N) content per leaf dry mass mg g-1 17.359583
13 158 Acer rubrum Leaf nitrogen (N) content per leaf dry mass mg N g-1 19.623538
14 181 Acer rubrum Absorptive fine root carbon/nitrogen (C/N) ratio Ratio 23.255814
In [31]:
Try[Try.trait=="Leaf nitrogen (N) content per leaf dry mass"]

print("Example:Leaf Nitrogen")
Try.head(n=20)
Example:Leaf Nitrogen
Out[31]:
Unnamed: 0 sp trait unit value
0 136 Acer pensilvanicum Leaf nitrogen (N) content per leaf dry mass mg g-1 24.200000
1 22 Acer pensylvanicum Leaf nitrogen (N) content per leaf dry mass % 2.412114
2 70 Acer pensylvanicum Leaf lignin content per leaf dry mass % dry mass 10.102437
3 177 Acer pensylvanicum Leaf nitrogen (N) content per leaf dry mass mg/g dry mass 26.258418
4 1 Acer rubrum Coarse woody debris (CWD) lignin content per C... NaN 77.750000
5 9 Acer rubrum Leaf nitrogen (N) content per leaf dry mass NaN NaN
6 16 Acer rubrum Branch coarse woody debris (CWD) nitrogen (N) ... % 0.264200
7 23 Acer rubrum Leaf nitrogen (N) content per leaf dry mass % 1.752372
8 59 Acer rubrum Stem (wood) carbon (C) content per stem dry mass % 48.640000
9 73 Acer rubrum Crown (canopy) width cm 2.812360
10 98 Acer rubrum Leaf carbon/nitrogen (C/N) ratio kg C per kg N 24.700000
11 111 Acer rubrum Leaf carbon (C) content per leaf area kg C per m2 39.900000
12 137 Acer rubrum Leaf nitrogen (N) content per leaf dry mass mg g-1 17.359583
13 158 Acer rubrum Leaf nitrogen (N) content per leaf dry mass mg N g-1 19.623538
14 181 Acer rubrum Absorptive fine root carbon/nitrogen (C/N) ratio Ratio 23.255814
15 2 Acer saccharum Coarse woody debris (CWD) lignin content per C... NaN 75.000000
16 10 Acer saccharum Leaf nitrogen (N) content per leaf dry mass NaN NaN
17 17 Acer saccharum Branch coarse woody debris (CWD) nitrogen (N) ... % 0.920000
18 24 Acer saccharum Leaf nitrogen (N) content per leaf dry mass % 1.819828
19 60 Acer saccharum Stem (wood) carbon (C) content per stem dry mass % 49.320000

These values were used to update values for N:C and lignin, These records can be found at in Cleanedup NECN documentation.csv. Values were also updated using records: primarily from Davis 2009

In [33]:
PATH = "E:/FS_meeting/FiguresMisc/"
Image(filename = PATH + "Davis2009.PNG", width=1000, height=1000)
Out[33]:

Values that were still missing were assesed as similar through genus then through genus then family and in a couple cases case order. These decesions can be found in the document:

GGDmin,GDDmax,Frost,D3,FRT

These values were taken from the original linkages manual which can be found at https://daac.ornl.gov/daacdata/model_archive/LINKAGES/comp/ORNL_TM-9519.pdf. Other values were taken from existing LANDIS papers, these can be found in the NECN folder. Species that could not be found in either way were adapted from qualitative assessment of range in comparison to know values for other species.

Functional Groups

In [36]:
importr("MASS", lib_loc= Where_is_R)
importr("meanShiftR", lib_loc= Where_is_R)
importr("ggfortify", lib_loc= Where_is_R)
importr("ggplot2", lib_loc= Where_is_R)
Out[36]:
rpy2.robjects.packages.Package as a <module 'ggplot2'>
In [37]:
%%R
library(meanShiftR)
library(ggfortify)
library(ggplot2)
#Set up the work drive
w_dir<-"E:/FS_meeting/Functional_Groups/"
setwd(w_dir)

To assess the functional groups I looked at the Ty Wilson maps for tree range and coupled them elevation levels, and 30yr normals for percipitation, temperature(min,mean,max), and vpd. These then used those to find the bounds of their range (maximum and minimum's of temp percip elevation, vpd). This was then used with a mean-shift algorthim to find the the closest groups based on mulitple criteria. Combinations were recombined to find groups. These groups were then plotted alongside two PCAs to display the grouping

In [38]:
%%R
read<-read.csv(paste(w_dir,"output.csv",sep=""),stringsAsFactors=FALSE)
Sampledf<-read[,-1]
Variableset<-c("sp","Type","Maxtemp","Mintemp","Meantemp","MinPPT","Minelevation","Maxelevation","MinVPD")
Variableset2<-c("sp","Type","Meantemp","MinPPT","MinVPD","Maxelevation")#Noelevation

Sampledf<-Sampledf[,Variableset2]

This is the conifer work

In [39]:
%%R
##Setting up the data for the bms.clustering
Clade="Conifer"
OneClade<-Sampledf[Sampledf[,2]==Clade,]
tree.labels<-OneClade[,c("sp")]
OneCladeproc<-((OneClade[c(-1,-2)]))
Height<-dim(OneCladeproc)[1]
Width<-dim(OneCladeproc)[2]

t<-as.numeric(data.matrix(OneCladeproc))


#seeds.data<-matrix(t,Width,Height)
row1<-t[c(1:Height)]
row1<- (row1-mean(row1))/sd(row1)
row2<-t[c((1+1*Height):(2*Height))]
row2<- (row2-mean(row2))/sd(row2)
row3<-t[c((1+2*Height):(3*Height))]
row3<- (row3-mean(row3))/sd(row3)
row4<-t[c((1+3*Height):(4*Height))]
row4<- (row4-mean(row4))/sd(row4)

seeds_data<-rbind(row1,row2,row3,row4)
seeds_data<-t(seeds_data)
names<-colnames(Sampledf)[c(-1,-2)]
In [44]:
%%R
b=.8

b=c(rep(b,ncol(seeds_data)))
iter<-10000
MS_Conif<- meanShift(
  seeds_data,
  seeds_data,
  bandwidth=b,
  alpha=0,
  iterations = iter
)
MS_Conif_assign <- MS_Conif$assignment
(print(MS_Conif_assign))
# value
MS_Conif_value <- MS_Conif$value
(print(MS_Conif_value))
##What is this
MSOUTPUT<-rbind(Variableset2[c(-1,-2)],MS_Conif_value)

par( mfrow=c( 1,2) )
for(i in 1:(length(names)-1)){
  #jpeg(paste(Clade,i,".jpg",sep=""))
  plot( seeds_data[,(i)], seeds_data[,(i+1)], col=MS_Conif$assignment,
        xlab=names[(i)], ylab=names[(i+1)], main=paste(Clade,ncol(seeds_data),"Groups",sep=""),
        cex=0.65, pch=16 )
  #dev.off()  
}
ms.labels6 <- MS_Conif$assignment
print( ms.labels6 )
ms.modes6 <- MS_Conif$value
print(ms.modes6)
B<-cbind(OneClade,ms.labels6)

seeds.plot<-data.frame((seeds_data))
colnames(seeds.plot)<-names
seeds.plot<-as.data.frame(cbind(seeds.plot,B$sp,as.character(B$ms.labels6)))
colnames(seeds.plot)[6]<-"MS"
colnames(seeds.plot)[5]<-"sp"
seeds.plot$MS<-as.factor(seeds.plot$MS)
seeds.plot$sp<-as.factor(seeds.plot$sp)

###Preform and plot PCA
PCA1<-prcomp(seeds.plot[,1:4])
#jpeg(paste(Clade,"PCA.jpg",sep=""))
autoplot(PCA1,data=seeds.plot,colour='MS',
         label=TRUE,loadings=TRUE,
         loading.label=TRUE,loadings.label.size=10,
         frame=TRUE
          )
#dev.off()

Hardwoods

In [45]:
%%R


read<-read.csv(paste(w_dir,"output.csv",sep=""),stringsAsFactors=FALSE)
Sampledf<-read[,-1]
Variableset<-c("sp","Type","Maxtemp","Mintemp","Meantemp","MinPPT","Minelevation","Maxelevation","MinVPD")
Variableset2<-c("sp","Type","Meantemp","MinPPT","MinVPD","Maxelevation")#Noelevation

Sampledf<-Sampledf[,Variableset2]

###Conifers#
Clade="Hardwoods"
OneClade<-Sampledf[Sampledf[,2]==Clade,]
tree.labels<-OneClade[,c("sp")]
OneCladeproc<-((OneClade[c(-1,-2)]))
Height<-dim(OneCladeproc)[1]
Width<-dim(OneCladeproc)[2]

t<-as.numeric(data.matrix(OneCladeproc))
print(t)

#seeds.data<-matrix(t,Width,Height)
row1<-t[c(1:Height)]
row1<- (row1-mean(row1))/sd(row1)
row2<-t[c((1+1*Height):(2*Height))]
row2<- (row2-mean(row2))/sd(row2)
row3<-t[c((1+2*Height):(3*Height))]
row3<- (row3-mean(row3))/sd(row3)
row4<-t[c((1+3*Height):(4*Height))]
row4<- (row4-mean(row4))/sd(row4)

seeds.data<-rbind(row1,row2,row3,row4)
seeds_data<-t(seeds.data)
names<-colnames(Sampledf)[c(-1,-2)]

##Set to three groups

b=.440

b=c(rep(b,ncol(seeds_data)))
iter<-100000
MS_HWD<- meanShift(
  seeds_data,
  seeds_data,
  bandwidth=b,
  alpha=0,
  iterations = iter
)
MS_HWD_assign <- MS_HWD$assignment
(print(MS_HWD_assign))
# value
MS_HWD_value <- MS_HWD$value
(print(MS_HWD_value))
##What is this
MSOUTPUT<-rbind(Variableset2[c(-1,-2)],MS_HWD_value)

par( mfrow=c( 1,2) )
for(i in 1:(length(names)-1)){
  #jpeg(paste(Clade,i,".jpg",sep=""))
  plot( seeds_data[,(i)], seeds_data[,(i+1)], col=MS_HWD$assignment,
        xlab=names[(i)], ylab=names[(i+1)], main=paste(Clade,ncol(seeds_data),"Groups",sep=""),
        cex=0.65, pch=16 )
  #dev.off()  
}


ms.labels6 <- MS_HWD$assignment
print( ms.labels6 )
ms.modes6 <- MS_HWD$value
print(ms.modes6)
B<-cbind(OneClade,ms.labels6)
In [47]:
%%R
seeds.plot<-data.frame(seeds_data)
colnames(seeds.plot)<-names
seeds.plot<-as.data.frame(cbind(seeds.plot,B$sp,as.character(B$ms.labels6)))
colnames(seeds.plot)[6]<-"MS"
colnames(seeds.plot)[5]<-"sp"
seeds.plot$MS<-as.factor(seeds.plot$MS)
seeds.plot$sp<-as.factor(seeds.plot$sp)

###Preform and plot PCA
PCA1<-prcomp(seeds.plot[,1:4])
#jpeg(paste(Clade,"PCA.jpg",sep=""))
autoplot(PCA1,data=seeds.plot,colour='MS',
         label=TRUE,loadings=TRUE,
         loading.label=TRUE,loadings.label.size=10,
         frame=TRUE
)
#dev.off()
print(PCA1)
PCAprint<-PCA1$rotation
EigenValues<-PCA1$sdev
PCAprint1<-rbind(EigenValues,PCAprint)
PCAprint1
#write.csv(PCAprint1,paste(Clade,"PCA1.csv",sep=""))

#write.csv(seeds.plot,paste(Clade,"FunctionalGroups.csv",sep=""))
#write.csv(MSOUTPUT,paste(Clade,"MS_Output.csv",sep=""))

I chose min ppt, min VPD, mean temp and mean elevation and ran a mean shift algorithm at different band widths until I found 3-4 groups. For the Conifers three groups fit well. For hardwoods there were a few species that fit no group. Because of how mean shift works expanding groups to include them collapsed too much of the difference between groups. In essence you can have 8 groups or 1 groups. I assigned these to the group closest too them, ending with 4 groups. List groups here.

Biomass

In this script I am looking to gain parameters around which to start parameterizing the biomass and growth curves in the NECN file of LANDIS-II for the Southern Apps Project. Max biomass is a hypotheical maximum that a plot could hold if not in competetion. To find that I am going to use known values of biomass and project them beyond a likly maximum. This file is a continuation of the a sorter that went through each FIA plot for the states of North Carolina, Tennessee, South Carolina, and Georgia, and calculated the Basal area for the plot and each species within it.

In [5]:
%%R
#Set up the work drive
w_dir<-"E:/Species Imput Files_1_3/Max_AGB/"
setwd(w_dir)

This is a look up table for the species in the study area. Here I am just grabbing out species names

In [6]:
%%R
SPLUT<-read.csv("SpeciesLUT.csv")
uniquesp<-SPLUT[,2]
print(length(SPLUT[,1]))
Important_Species<-c(129,621,832,316,806,833,132,711,802,261)

In this following loop we will

  • Calculate the percent of total biomass
  • Calculate the ratio of species biomass to total biomass
  • From the upper half of this ratio we will preform a linear regression
  • Where this regression is at a 150% occupany of a species, we will set as Max AGB
  • Review graphs of this process and save the results
In [22]:
%%R
#par(mfrow=(2,20))
df<-NULL

for(species in uniquesp){
  #print(species)
 Speciesname<-as.character(SPLUT[,13][SPLUT[,2]==species])

  read<-read.csv(paste(w_dir,species,"AGB.csv",sep=""))
  read<-as.data.frame(read[complete.cases(read),])
  #uniquepercent<-unique(read$perecentofAGB)
 # for(t in uniquepercent){
  #  agebreak<-read[read$perecentofAGB==t,]
   # agebreak<-agebreak[agebreak$`AboveGround Carbon`> quantile(agebreak$`AboveGround Carbon`,probs=.80),]
  #  out<-rbind(agebreak,out)
  #}
  

  read$ration<-(read$sumofsp/read$perecentofAGB)
  read$ration[is.na(read$ration)]<-0.0
  readquant<-read[read$ration > quantile(read$ration,probs=.75),]
  readquant$adjsum<-readquant$sumofsp*0.112085
  reg1<-lm(readquant$adjsum~readquant$perecentofAGB)
  #max<-max(read$sumofsp[!is.na(read$sumofsp)])
  maxplot<-max(readquant$sumofsp*0.112085)
  #SteepestMax<-maxplot$sumofsp/maxplot$perecentofAGB
  
  intercept<-as.numeric(reg1$coefficients[1])

  slope<-as.numeric(reg1$coefficients[2])
  likelymaximum<-((1.00)*slope)+intercept
  #print(likelymaximum)#to## g per m2
  jpeg(filename=(paste(species,"_svp.jpeg",sep="")))
  plot(readquant$adjsum~readquant$perecentofAGB, main=paste(Speciesname,",","n=",length(readquant$sumofsp),",","Max AGB",round(likelymaximum,2))
       ,xlim=c(0,2.0),ylim=c(0,likelymaximum*1.5),xlab="Percent of Stand",ylab="Above ground biomass g/m2",col="Green",pch = 19,cex.main=2.0,cex.lab=1.4,cex.axis=1.2,cex=1.5,font=2)
  abline(reg1)
  abline(v=1.0,col="red")
  abline(h=likelymaximum,col="red")
  dev.off()
    
  intercept<-as.numeric(reg1$coefficients[1])
  slope<-as.numeric(reg1$coefficients[2])

  #SteepestMax*.112085
A<-cbind(species,likelymaximum,maxplot,slope,intercept)
df<-rbind(A,df)
}

colnames(df)<-c("Species","Linear Regression(gm-2) at 125% ","Maximum Measured(gm-2)","Slope","Intercept")
write.csv(w_dir,"AGBlist.csv")
In [23]:
from matplotlib.pyplot import figure, imshow, axis
from matplotlib.image import imread

mypath='E:/Species Imput Files_1_3/Max_AGB'

hSize = 30
wSize = 50
col = 3

Important_Species=[129,621,832,316,806,833,132,711,802,261]

Growthcurveplot=["129_svp.jpeg","621_svp.jpeg","832_svp.jpeg","832_svp.jpeg","316_svp.jpeg","806_svp.jpeg","833_svp.jpeg",
                 "132_svp.jpeg","711_svp.jpeg","802_svp.jpeg","261_svp.jpeg",]

showImagesMatrix(Growthcurveplot, hSize = 62, wSize = 20,col=2)

Growth curves

Now that we have seen what some values for the maximum amount of Biomass we want to get a better idea for what the growth curves should look like so that we can calibarte the growth over time. To do this we use a "" associaation to look at the sight index and height to relate to age. Then we will look at the upper 20% of sites a fit a logarithmic regression to simulate the growht patterns of trees. This will be used as a basis for LANDIS-II Trials to tune the rate of growth.

Importing libraries from R

In [ ]:
importr("plyr", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
#importr("dplyr", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
importr("sp", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
#importr("magrittr", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")

Setting the Working Drive

In [ ]:
%%R
library(plyr)
#library(dplyr)
library(sp)
library(raster); options(scipen=999) # scipen removes scientific notation from FIA data

#library(rgdal)

wdir<-"C:/Users/zjrobbin/Desktop/Growth_cruves/"
setwd(wdir)
In [ ]:
%%R
TPA.factor <- 6.018046 #trees/acre adjust (TPA in FIA)
conv.factor <- .112085 #Convert from lbs/acre to g/m-2
In [ ]:
%%R
###The Output from IC 2 Get Age
AgeGraphs<-read.csv(paste(wdir,"All_Sapps_FIA_AGE_TREE_HT.csv",sep=""))
print(head(AgeGraphs))
Plts<-unique(AgeGraphs$PLT_CN)
unique(AgeGraphs$calc_age_rounded)
Plts<-unique(AgeGraphs$PLT_CN)




colnames(AgeGraphs)
AgeGraphs$AdjCarbon<-AgeGraphs$CARBON_AG*AgeGraphs$TPAGROW_UNADJ*conv.factor
output<-NULL
In [ ]:
%%R



for(i in 1:length(Plts)){
One<-as.data.frame(AgeGraphs[AgeGraphs$PLT_CN==Plts[i],])
plt<-Plts[i]
IC<-aggregate(round(One$AdjCarbon),by=list(SPCD=One$SPCD,AGE=as.numeric(One$calc_age_rounded)),FUN=sum)
IC$PLT_CN<-plt
colnames(IC)<-c("SPCD","AGE","AboveGround Carbon","PLT_CN")
output<-rbind(IC,output)
}

print(output)
In [ ]:
%%R
print(output)
write.csv(output,paste(wdir,"FIA_Height_Carbon_Aggregate.csv",sep=""))
In [ ]:
%%R

output<-read.csv(paste(wdir,"FIA_Height_Carbon_Aggregate.csv",sep=""))
print(head(output))
SPLUT<-read.csv("SpeciesLUT.csv")
uniquesp<-SPLUT[,2]
uni.sp<-unique(output$SPCD)
output<-output[output$AboveGround.Carbon>0,]
#print(colnames(SPLUT))
print(colnames(output))

This is the label for the for loop that lies below as the ...

In [ ]:
%%R

for(i in uni.sp){
  IC<-output[output$SPCD==i,]
  #print("Species")
  #print(i)
  Speciesname<-as.character(SPLUT[,13][SPLUT[,2]==i])
  out<-NULL
  if(length(IC$AGE)>100){
  uniqueage<-unique(IC$AGE)
  
  for(t in uniqueage){
  agebreak<-IC[IC$AGE==t,]
  agebreak<-agebreak[agebreak$AboveGround.Carbon> quantile(agebreak$AboveGround.Carbon,probs=.80),]
  out<-rbind(agebreak,out)
  }
  }
  if(!is.null(out)){
  IC<-out
  }
  #print(length(IC$AGE))
  x<-IC$AGE
  y<-IC$AboveGround.Carbon
  
  filename<-paste(wdir,"Output/",i,"_GY_Curve.jpeg",sep="")
  #print(filename)
  jpeg(file = filename, bg = "transparent")
  #plot(1:10)
  #rect(1, 5, 3, 7, col = "white")
  title<-paste(Speciesname, i, "n=",length(IC$AGE),"plots")
  plot(y~x,main=title,xlim=c(0,200),ylim=c(0,4000),xlab="Age",ylab="Carbon g/m2",pch = 16,col="gray",cex.main=2.0,cex.lab=1.5,cex.axis=1.2,font=2)
  fit<-lm(y~log(x))
  summary(fit)
  slope<-as.numeric(fit$coefficients[2])
  lm(formula=y~log(x))
  n=500
  set.seed(10)
  x=seq(from=0,to=200,length.out=201)
  y=predict(fit,newdata=list(x=seq(from=0,to=200,length.out=201)),
            interval="confidence")
  g_m2_year<-(slope/x)*5
  write<-(cbind(x,y,g_m2_year))
 
  matlines(x,y,lwd=2)
  dev.off()
  x<-IC$AGE
  y<-IC$AboveGround.Carbon
    ##Just plot important Species

    
  write.csv(write,paste(wdir,"Output/",i,"GR_.csv",sep=""))
}

Here are graphs the 9 most prevelant species on the landscape. For each functional group we will use the two most prevalent species and find parameters that best align with the growth curve.

In [49]:
from matplotlib.pyplot import figure, imshow, axis
from matplotlib.image import imread

mypath='E:/Growth_cruves/Output'

hSize = 30
wSize = 50
col = 3

Important_Species=[129,621,832,316,806,833,132,711,802,261]

Growthcurveplot=["129_GY_Curve.jpeg","621_GY_Curve.jpeg","832_GY_Curve.jpeg","832_GY_Curve.jpeg","316_GY_Curve.jpeg","806_GY_Curve.jpeg","833_GY_Curve.jpeg","132_GY_Curve.jpeg",
                "711_GY_Curve.jpeg","802_GY_Curve.jpeg","261_GY_Curve.jpeg",]

showImagesMatrix(Growthcurveplot, hSize = 62, wSize = 20,col=2)

Aboitic Model Single Cells: Probably Will have to redo

In [51]:
#Carbon
PATH = "E:/FS_meeting/FiguresMisc/"
Image(filename = PATH + "Carbon.jpg", width=1000, height=1000)
Out[51]:
In [52]:
#Nitrogen
PATH = "E:/FS_meeting/FiguresMisc/"
Image(filename = PATH + "Nitrogen.jpg", width=1000, height=1000)
Out[52]:
In [53]:
##Other Ecosystem stuff
mypath='E:/FS_meeting/FiguresMisc'

Growthcurveplot=["LAI.jpg","Eco.jpg"]

showImagesMatrix(Growthcurveplot, hSize = 100, wSize = 70,col=2)

Functional Groups Parameterizing: Ongoing

In [55]:
#importr("readtext", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
importr("scales", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
importr("RColorBrewer", lib_loc= "C:/Users/zjrobbin/Documents/R/win-library/3.4")
Out[55]:
rpy2.robjects.packages.Package as a <module 'RColorBrewer'>
In [56]:
%%R
#library(readtext)
library(scales)
library(RColorBrewer)
LandisDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/LANDIS_Sapps_Singlecell_Trial_SS/"
setwd(LandisDrive)##Wdir has to be landisDrive
writeDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/SpeciesFiles/"
DataDrive<-"C:/Users/zjrobbin/Desktop/Growth_cruves/Output/"
display.brewer.all()
Error in setwd(LandisDrive) : cannot change working directory
---------------------------------------------------------------------------
RRuntimeError                             Traceback (most recent call last)
~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\ipython\rmagic.py in eval(self, code)
    253             # Need the newline in case the last line in code is a comment
--> 254             value, visible = ro.r("withVisible({%s\n})" % code)
    255         except (ri.RRuntimeError, ValueError) as exception:

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\robjects\__init__.py in __call__(self, string)
    351         p = _rparse(text=StrSexpVector((string,)))
--> 352         res = self.eval(p)
    353         return conversion.ri2py(res)

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\robjects\functions.py in __call__(self, *args, **kwargs)
    177                 kwargs[r_k] = v
--> 178         return super(SignatureTranslatedFunction, self).__call__(*args, **kwargs)
    179 

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\robjects\functions.py in __call__(self, *args, **kwargs)
    105             new_kwargs[k] = conversion.py2ri(v)
--> 106         res = super(Function, self).__call__(*new_args, **new_kwargs)
    107         res = conversion.ri2ro(res)

RRuntimeError: Error in setwd(LandisDrive) : cannot change working directory


During handling of the above exception, another exception occurred:

RInterpreterError                         Traceback (most recent call last)
~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\ipython\rmagic.py in R(self, line, cell, local_ns)
    703             else:
--> 704                 text_result, result, visible = self.eval(code)
    705                 text_output += text_result

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\ipython\rmagic.py in eval(self, code)
    256             warning_or_other_msg = self.flush() # otherwise next return seems to have copy of error
--> 257             raise RInterpreterError(code, str(exception), warning_or_other_msg)
    258         text_output = self.flush()

RInterpreterError: Failed to parse and evaluate line '#library(readtext)\nlibrary(scales)\nlibrary(RColorBrewer)\nLandisDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/LANDIS_Sapps_Singlecell_Trial_SS/"\nsetwd(LandisDrive)##Wdir has to be landisDrive\nwriteDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/SpeciesFiles/"\nDataDrive<-"C:/Users/zjrobbin/Desktop/Growth_cruves/Output/"\ndisplay.brewer.all()\n'.
R error message: 'Error in setwd(LandisDrive) : cannot change working directory'

During handling of the above exception, another exception occurred:

PermissionError                           Traceback (most recent call last)
<ipython-input-56-7a67d1c5c5df> in <module>
----> 1 get_ipython().run_cell_magic('R', '', '#library(readtext)\nlibrary(scales)\nlibrary(RColorBrewer)\nLandisDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/LANDIS_Sapps_Singlecell_Trial_SS/"\nsetwd(LandisDrive)##Wdir has to be landisDrive\nwriteDrive<-"C:/Users/zjrobbin/Desktop/Sapps_SC/SpeciesFiles/"\nDataDrive<-"C:/Users/zjrobbin/Desktop/Growth_cruves/Output/"\ndisplay.brewer.all()\n')

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
   2319             magic_arg_s = self.var_expand(line, stack_depth)
   2320             with self.builtin_trap:
-> 2321                 result = fn(magic_arg_s, cell)
   2322             return result
   2323 

<decorator-gen-130> in R(self, line, cell, local_ns)

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188 
    189         if callable(arg):

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\rpy2-2.9.4-py3.6-win-amd64.egg\rpy2\ipython\rmagic.py in R(self, line, cell, local_ns)
    717                 print(e.err)
    718             if tmpd:
--> 719                 rmtree(tmpd)
    720             return
    721         finally:

~\Anaconda3\envs\Python and R 1_29\lib\shutil.py in rmtree(path, ignore_errors, onerror)
    492             os.close(fd)
    493     else:
--> 494         return _rmtree_unsafe(path, onerror)
    495 
    496 # Allow introspection of whether or not the hardening against symlink

~\Anaconda3\envs\Python and R 1_29\lib\shutil.py in _rmtree_unsafe(path, onerror)
    387                 os.unlink(fullname)
    388             except OSError:
--> 389                 onerror(os.unlink, fullname, sys.exc_info())
    390     try:
    391         os.rmdir(path)

~\Anaconda3\envs\Python and R 1_29\lib\shutil.py in _rmtree_unsafe(path, onerror)
    385         else:
    386             try:
--> 387                 os.unlink(fullname)
    388             except OSError:
    389                 onerror(os.unlink, fullname, sys.exc_info())

PermissionError: [WinError 32] The process cannot access the file because it is being used by another process: 'C:\\Users\\zacha\\AppData\\Local\\Temp\\tmpz0kwuwf0\\Rplots001.png'

Using the the growth curves and the values of previous landis functional groups as a estimate we will paramterize the functional groups. Here is an example of the proccesing Quercus Prinus(sp=832) in the Northern Hardwood Group. Using the minimum positive value from the growth curve and its agb value, I grow one cohort of that value with the establishment turned to zero. First we ajdust the LAI factors.

In [57]:
%%R
#Paramaters
spp<-"129"

Functional Group LAI parameters: Affect AGB considerably and need to parameterized first.

  • BTOLAI: Affects very little of LAI Total.
  • KLAI: Affects rate of curve.
  • MAXLAI: Affects scaler of the curve.
  • </ul> These are compared with functional group parameters from Hardiman et al., 2013

In [ ]:
%%R

Log<-"NECN-calibrate-log.csv"
Cal<-read.csv(paste(LandisDrive,Log,sep=""),row.names = NULL)

Cal$redux<-as.numeric(Cal$row.names)+10
Lastcal<-Cal###Remove after demo
plot(Cal$LAI~Cal$redux,main="LAI Total",xlab="Years", ylab="Total LAI",xlim=c(0,100),ylim=c(0,8),type="p",pch = 19,cex.main=2.0,cex.lab=1.35,cex.axis=1.2,cex=1.5,font=15,font.lab=15)

plot(Lastcal$LAI~Lastcal$redux,main="Previous Run",xlab="Years", ylab="Total LAI",xlim=c(0,100),ylim=c(0,8),,type="p",pch = 19,cex.main=2.0,cex.lab=1.45,cex.axis=1.2,cex=1.5,font=15,font.lab=15)
In [58]:
PATH = "C:/Users/zjrobbin/Desktop/Sapps_SC/PicsofTrees/"
Image(filename = PATH + "Oak_Hickory.png", width=500, height=500)
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-58-59d5de66f720> in <module>
      1 PATH = "C:/Users/zjrobbin/Desktop/Sapps_SC/PicsofTrees/"
----> 2 Image(filename = PATH + "Oak_Hickory.png", width=500, height=500)

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\display.py in __init__(self, data, url, filename, format, embed, width, height, retina, unconfined, metadata)
   1161         self.unconfined = unconfined
   1162         super(Image, self).__init__(data=data, url=url, filename=filename, 
-> 1163                 metadata=metadata)
   1164 
   1165         if self.width is None and self.metadata.get('width', {}):

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\display.py in __init__(self, data, url, filename, metadata)
    607             self.metadata = {}
    608 
--> 609         self.reload()
    610         self._check_data()
    611 

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\display.py in reload(self)
   1192         """Reload the raw data from file or URL."""
   1193         if self.embed:
-> 1194             super(Image,self).reload()
   1195             if self.retina:
   1196                 self._retina_shape()

~\Anaconda3\envs\Python and R 1_29\lib\site-packages\IPython\core\display.py in reload(self)
    632         """Reload the raw data from file or URL."""
    633         if self.filename is not None:
--> 634             with open(self.filename, self._read_flags) as f:
    635                 self.data = f.read()
    636         elif self.url is not None:

FileNotFoundError: [Errno 2] No such file or directory: 'C:/Users/zjrobbin/Desktop/Sapps_SC/PicsofTrees/Oak_Hickory.png'

Citation: Hardiman, B. S., Gough, C. M., Halperin, A., Hofmeister, K. L., Nave, L. E., Bohrer, G., & Curtis, P. S. (2013). Maintaining high rates of carbon storage in old forests: a mechanism linking canopy structure to forest function. Forest Ecology and Management, 298, 111-119.

Once these are in the range of correct we can work on tuning the aboveground biomass. Name<-"NHW" Index<-"4" PPDF1<-"25.0" PPDF2<-"40.0" PPDF3<-"2.05" ###beginning arc higher is lower .5-.05 PPDF4<-"5.0"###long arc Higher is lower whole number FCFRAC<-"0.48"

  • PPDF1: The minimum temperature of the range
  • PPDF2: The maximum temperature of the range
  • PPDF3: The earlier trajectory of the curve of AGB
  • PPDF4: The later trajectory of the curve of AGB
  • </ul> This is using the growth curves we put together earlier.

    Using each run and an ealier run we adjust the parameters to best match the aboveground AGB

In [ ]:
%%R
RunName<-"QuerPrin4"
Log<-"NECN-succession-log.csv"


FileWrite<-paste(writeDrive,RunName,sep="")
#dir.create(FileWrite)
print("set Age at carbon")

NECNSuccession<-read.csv(paste(LandisDrive,Log,sep=""))
GrowthCurves<-read.csv(paste(DataDrive,spp,"GR_.csv",sep=""))




        ##
NECNSuccession<-read.csv(paste(LandisDrive,Log,sep=""))
GrowthCurves<-read.csv(paste(DataDrive,spp,"GR_.csv",sep=""))
Growthxy<-GrowthCurves[,c(2,3)]
Growthxy$fit<-Growthxy$fit*2
NECNxy<-NECNSuccession[,c(1,7)]
Start<-Growthxy[Growthxy$fit>0,][1,]
Agestart<-Start[1,1]
NECNxy$Time<-NECNSuccession$Time+Agestart
#print(Start)
#par(mfrow=c(1,1))
#dev.new(width=5, height=4)
plot(NECNxy,col=alpha("darkorange"),xlim=c(10,120),ylim=c(0,10000),type="p",main="This Run",pch = 19,cex.main=2.0,cex.lab=1.4,cex.axis=1.2,cex=1.5,font=2)
points(Growthxy,pch = 19,type="p",col=alpha("darkgreen",0.8),cex=1.5)
legend("topleft",legend=c("Landis","FIA"),col=c("darkorange","darkgreen"),cex=1.0,pt.cex=1.0,pch=c(19,19))

filename<-paste(paste("Spp",spp,".png",sep=""))
#print(filename)

plot(PreviousRunNECN,col=alpha("darkorange"),xlim=c(10,120),ylim=c(0,10000),type="p",main="Previous Run",pch = 19,cex.main=2.0,cex.lab=1.4,cex.axis=1.2,cex=1.5,font=2)
points(PreviousRunGrowth,pch = 19,type="p",col=alpha("darkgreen",0.8),cex=1.5)
legend("topleft",legend=c("Landis","FIA"),col=c("darkorange","darkgreen"),cex=1.0,pt.cex=1.0,pch=c(19,19))
In [ ]:
%%R
PreviousRunGrowth<-Growthxy
PreviousRunNECN<-NECNxy

jpeg(file =paste(writeDrive,"/",RunName,"/",RunName,".jpeg",sep=""), bg = "transparent")
plot(NECNxy,col="red",xlim=c(16,120),ylim=c(0,10000))
points(Growthxy)
legend("topleft",legend=c("Landis","FIA"),col=c("red","black"),cex=1.0,pt.cex=1.0,pch=c(1,1))
dev.off()

object<-readtext(paste(LandisDrive,"Sapps_NECN_Values_SingleCell.txt",sep=""))

fileConn<-file(paste(writeDrive,RunName,"/",RunName,"LandscapeText.txt",sep=""))
writeLines(object[,2], fileConn)
close(fileConn)

write.csv(NECNSuccession,paste(writeDrive,RunName,"/",RunName,"sucesssion.csv",sep=""))

We can use the same data as before to compare to the NPP from our LANDIS run

In [ ]:
%%R
NEEC<-cbind(NECNSuccession[,(1)],(NECNSuccession[,(8)]+NECNSuccession[,(9)])*1/90.7185)
plot(NEEC,xlim=c(16,120),ylim=c(0,10),type="p",pch = 19,cex.main=2.0,cex.lab=1.45,cex.axis=1.2,cex=1.5,font=15,font.lab=15,,main="NPP t C Ha-1 yr-1",xlab="Years", ylab="tC HA-1 yr-1")
In [ ]:
PATH = "C:/Users/zjrobbin/Desktop/Sapps_SC/PicsofTrees/"
Image(filename = PATH + "Oak_Hickory.png", width=500, height=500)

Take away: